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Abstract 

The IceCube Neutrino Observatory, approximately 1 km^ in size, is now complete with 86 
strings deployed in the Antarctic ice. IceCube detects the Cherenkov radiation emitted by 
charged particles passing through or created in the ice. To realize the full potential of the de- 
tector, the properties of hght propagation in the ice in and around the detector must be well 
understood. This report presents a new method of fitting the model of light propagation in the 

2 



ice to a data set of in-situ light source events collected with IceCube. The resulting set of derived 
parameters, namely the measured values of scattering and absorption coefficients vs. depth, is 
presented and a comparison of IceCube data with simulations based on the new model is shown. 



1. Introduction 

IceCube is a cubic-kilometer-scale high-energy neutrino observatory built at the geographic 
South Pole |7] (see Fig. [T]). A primary goal of IceCube is to elucidate the mechanisms for produc- 
tion of high-energy cosmic rays by detecting high-energy neutrinos from astrophysical sources. 
IceCube uses the 2.8 km thick glacial ice sheet as a medium for producing Cherenkov light 
emitted by charged particles created when neutrinos interact in the ice or nearby rock. Neutrino 
interactions can create high-energy muons, electrons or tau leptons, which must be distinguished 
from a background of downgoing atmospheric muons based on the pattern of emitted Cherenkov 
light. This light is detected by an embedded array of 5160 optical sensors (digital optical mod- 
ules, or DOMs for short), 4680 of which are deployed at depths of 1450 - 2450 m and spaced 
17 m apart along 78 vertical cables (strings). The strings are arranged in a triangular lattice with 
a horizontal spacing of approximately 125 m. The remaining 480 sensors are deployed in a more 
compact geometry forming the center of the DeepCore array [2]. The IceCube optical sensors 
are remotely-controlled autonomous detection units which digitize the data. They include light- 
emitting diodes (LEDs) which may be used as artificial in-situ light sources. Also shown in Fig. 
\T\ is the location of the AMANDA-II neutrino telescope. AMANDA-II was the precursor for 
IceCube and was composed of 677 optical sensors organized along 19 strings, with most of the 
sensors located at depths of 1500 to 2000 m. It operated as a part of the IceCube observatory 
until it was decommissioned in May 2009. 

Cherenkov photons are emitted with a characteristic wavelength dependence of in the 
wavelength range of 300-600 nm, which includes the relevant sensitivity region of the photo- 
sensors. Photons are emitted in a cone around the direction of particle motion with an opening 
angle, determined by the speed of the particle and refractive index of the ice Ist], of about 41° for 
relativistic particles. As the photons propagate from the point of emission to the receiving sen- 
sor, they are affected by absorption and scattering in the ice. These propagation effects must be 
considered for both simulation and reconstruction of IceCube data and thus need to be carefully 
modeled. The important parameters to describe photon propagation in a transparent medium are: 
the average distance to absorption, the average distance between successive scatters of photons, 
and the angular distribution of the new direction of a photon at each given scattering point. This 
work presents a new, global-fit approach which achieves an improved description of experimental 
data. 

To determine the ice parameters, dedicated measurements are performed with the IceCube 
detector Photons are emitted by the LEDs in DOMs and recorded by other DOMs, as sketched 
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Figure 1 : The IceCube Neutrino Observatory, final configuration. Also shown is the AMANDA array, pre- 
cursor to IceCube, which ended operation in 2009. 




1000 2000 3000 4000 5000 
time from the flasher event [ ns ] 



Figure 2: Left (a): simplified schematics of the experimental setup: the flashing sensor on the left emits 
photons, which propagate through ice and are detected by a receiving sensor on the right. Right (b): example 
photon arrival time distributions at a sensor on one of the nearest strings (122 m away) and on one of the 
next-to-nearest strings (217 m away; histogram values are multiplied by a factor of 10 for clarity). Dashed 
lines show data and solid lines show simulation based on the model of this work (with best fit parameters). 
The goal of this work is to find the best-fit ice parameters that describe these distributions as observed in data 
simultaneously for all pairs of emitters and receivers. 



in Fig.|2^. The recorded data include the total charge (corresponding to the number of arriving 
photons) and photon arrival times, shown in Fig. |2]3. A data set that covers all detector depths 
was produced. A global fit of these data was performed, and the result is a set of scattering and 
absorption parameters that best describes the full data set. The AMANDA Collaboration used 
an analysis based on separate fits to data for individual pairs of emitters and receivers to 
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measure the optical properties of the ice. These fits used data taken at very low light levels, to 
avoid multi-photon pileup detector effects. 

The relevant detector instrumentation is described in Section |2] of this paper Section |3] in- 
troduces the data set. The parameterization for modeling the ice surrounding the detector is 
described in Section |4] while Section |5] discusses the simulation. The likelihood function used 
to compare data and simulation is discussed in Sections |6] and |7] and Section |8]explains how the 
search for the best solution was performed. Section |9]compares the result with an independent 
probe of the dust concentration in ice [5J. Finally, Section [10] discusses the uncertainties of the 
measurement, Section [TT]presents data-simulation comparisons, and Section [T2l summarizes the 
result. 

2. Instrumentation 

The data for this analysis were obtained in 2008 when IceCube consisted of 40 strings, each 
with 60 DOMs, as shown in Fig. [3] Each of the DOMs consists of a 10" photomultiplier tube 
(PMT) 10] facing downwards and several electronics boards enclosed in a glass pressure sphere 
yj]. The main board of the electronics includes two types of digitizers for recording PMT wave- 
forms as well as time stamping, control and communications (|7|]. The first 427 ns of each wave- 
form is digitized at 300 megasamples per second by fast ATWD chips (analog transient waveform 
digitizer, see [7J), and longer duration signals are recorded at 25 megasamples per second by the 
fast ADC (fast analog-to-digital converter, or fADC for short) chips. The system is capable of 
resolving charge of up to 300 photoelectrons per 25 ns with precision limited only by the prop- 
erties of the PMTs (i.e., 1 photoelectron is resolved with ~ 25% uncertainty in charge). Both the 
ATWD and fADC use 10 bits for amplitude digitization. However, the ATWD uses three parallel 
channels with different gains (with a factor of about 8 between) and has a finer time resolution 
than the fADC (roughly 3.3 vs. 25 ns bin width). The main board contains two ATWD chips on 
each DOM, ensuring that a waveform can be recorded with one chip while the other one is read 
out, thus reducing the sensor dead time. 

Each DOM includes 12 LEDs on a "flasher board" that produce pulsed light detectable by 
other DOMs located up to 0.5 km away. The primary purpose of the measurements with these 
flashers is calibration of the detector These calibration studies include determining the detector 
geometry, verifying the calibration of time offsets and the time resolution, verifying the linearity 
of photon intensity measurement, and extracting the optical properties of the detector ice (this 
paper). 

Depending on the intended application, flasher pulses can be programmed with rates from 
1.2 Hz to 610 Hz, durations of up to 70 ns, and LED currents up to 240 mA. The corresponding 
total output from each LED ranges from below 10^ to about lO'" photons. The programmed 
current pulse is applied to each individual LED through a high-speed MOSFET (metal-oxide- 
semiconductor field-effect transistor) driver with a series resistor. The voltage across the resistor 
is recorded by the DOM's waveform digitizer to precisely define the onset of each pulse. Figure|4] 
shows laboratory measurements of the optical output time profiles from short and wide pulses. 
Pulses exhibit a characteristic rise time of 3-4 ns and a small afterglow, decaying with a 12 ns 
time constant. The narrowest pulses achievable have a full width at half maximum (FWHM) of 
6 ns. 

The wavelength spectrum has been measured for the LED light exiting the glass pressure 
sphere and was found to be centered at 399 nm with a FWHM of 14 nm (see Fig. |5]). This 
wavelength was chosen to approximate the typical wavelength of detected Cherenkov photons 
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Figure 3: Left: Top view layout of IceCube in tlie 40-string configuration in 2008. String 63, for whicli tlie 
DOMs emitted flashing light in the study presented here, is shown in black. The nearest 6 strings are shown 
in brown. The dashed lines and numbers 2009 and 2010 in the left figure indicate the approximate location 
of the detector parts deployed during those years. Right: a typical DOM flasher event, DOM 46 on string 63 
flashing. The larger circles represent DOMs that recorded larger numbers of photons. The arrival time of the 
earhest photon in each DOM is indicated with color: eai'ly times are shown in red while late times trend to 
blue. 




Time (nsec) 

Figure 4: Flasher light output time profile for pulses of minimum and maximum width. The relative height of the 
short pulse has been scaled so the leading edges are comparable. This measurement was performed using a small PMT 
(Hamamatsu R1450) after optical attenuation of the pulses to facilitate counting of individual photons. 



(as discussed and shown in Fig. [8] below). To supplement data from the standard flashers, 16 
special DOMs were constructed and deployed with LEDs that emit at 340 nm, 370 nm, 450 nm, 
and 500 nm. Data from these special flashers were not used in the analysis of this paper but will 
be used in future analyses of wavelength-dependent effects. 

The 12 LEDs in each DOM are aimed in six different azimuth angles (with 60° spacing) and 
along two diff'erent zenith angles. After correcting for refraction at interfaces between air, glass 
and ice, the angular emission profiles peak along the horizontal direction for the 6 horizontal 
LEDs and 48° above the horizontal for the 6 tilted LEDs. The angular spread is reduced by the 
refraction and is modeled using a 2-D Gaussian profile with cr = 10° around each peak direction. 
During the DOM deployment and freeze-in within the glacial ice sheet, the azimuthal orientations 
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LED output spectrum at MB temperature -1 5°C 




Figure 5: Wavelength spectrum of light emitted for a DOM operating at a mainboard (MB) temperature of -15° C. The 
y-axis shows the average number of photons detected per pulse of the LED light. 



of the DOMs are not controlled and are initially unknown. The orientation of each DOM, and 
therefore the initial direction of emitted light from each LED, is determined to a precision of 
about 10° by flashing individual horizontal LEDs and studying the light arrival time at the six 
surrounding strings. Here one relies on direct light from an LED facing a target arriving sooner 
than scattered light from one facing away. 

3. Flasher data set 

The data set used in this paper includes at least 250 flashes from each DOM on string 63. 
DOMs were flashed at 1 .2 Hz in a sequence, using a 70 ns pulse width and maximum brightness. 
The six horizontal LEDs on each flasher board were operated simultaneously, creating a pattern 
of light with approximate azimuthal symmetry around the flasher string. Flash sequences for 
DOMs at dififerent depths were overlapping but were sufficiently displaced in time that pulses of 
observed light were unambiguously assigned to individual flashers. 

As seen in Fig.|6l there is a substantial variation among the charges collected in DOMs at ap- 
proximately the same depth as the emitter on the six surrounding strings. Some of the variation is 
due to relative differences in light yield between LEDs, and some is due to differences in distance 
to, and depth of, the six surrounding strings. Other reasons may include non-homogeneity of the 
ice. 

The pulses corresponding to the arriving photons were extracted from the digitized wave- 
forms and binned in 25 ns bins, from to 5000 ns from the start of the flasher pulse (extracted 
from the special-purpose ATWD channel of the flashing DOM). To reduce the contribution from 
saturated DOMs (most of which were near the flashing DOM on string 63) (§], and to minimize 
the effects of the systematic uncertainty in the simulated angular sensitivity model of a DOM, 
the photon data collected on string 63 were not used in the fit. A DOM becomes saturated when 
it is hit by so many photons that the charge in its digitized output is no longer proportional to the 
number of incident photons. 

The angular sensitivity model specifies a fraction of photons that are detected at a given angle 
with respect to the PMT axis. It accounts for the nominal DOM sensitivity measured in the lab, 
modified by the scattering in the column of re-frozen ice (see Fig. |7] and further discussion in 
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Figure 6: Charge collected by DOMs on the six nearest strings (121.8 - 126.6 m away, triangles) and six 
next-to-nearest strings (21 1.4- 217.9 m away, circles), observed when flashing at the same position on stiing 
63. 



section|5]l. Variations in the angular sensitivity model have a large impact on the simulated DOM 
response to the photons arriving along the PMT axis (straight into the PMT or into the back of a 
DOM), while the response to photons arriving from the sides of the PMT is much less affected. 

4. Six-parameter ice model 

This section overviews the ice parameterization introduced in f7], which in this paper is 
referred to as the six-parameter ice model. The ice is described by a table of depth-dependent 
parameters /7p(400) and adust(400) related to scattering and absorption at a wavelength of 400 nm, 
by the depth-dependent relative temperature 6t, and by the six global parameters (measured in 
||4|]): a, K, A, B, D, and E, which are described below. The thickness of the ice layers was 
somewhat arbitrarily chosen to be 10 m. The scattering and absorption coefficients of each ice 
layer are best interpreted as the average of their true values over the thickness of the ice layer 
The chosen thickness of 10 m is the same as the value chosen in |4] but smaller than the vertical 
DOM spacing of 17 m. Due to small depth offsets between the DOMs on different strings, we 
retain at least 1 receiving DOM per layer 

The geometrical scattering coefficient b determines the average distance between successive 
scatters (as l/b). It is often more convenient to quote the effective scattering coefficient, be - 
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b ■ (1 - {cos 6)), where is the deflection angle at each scatter. The absorption coefficient a 
determines the average distance traveled by a photon before it is absorbed (as I /a). 

The wavelength dependence of the scattering and absorption coefficients is given by the fol- 
lowing expressions (for wavelength A in nm). The power law dependence is predicted by theo- 
retical models of light scattering in dusty ice. The power law dependence on photon wavelength 
was verified in the AMANDA study, using light sources with several different frequencies 
The effective scattering coefficient, with the global fit parameter a, is 

tAA)^bJ400)-(^)\ 

The total absorption coefficient is the sum of two components, one due to dust and the other a 
temperature dependent component for pure ice |3| : 



aiA) = fldust(^) + Ae-*/'' ■ (1 + 0.01 ■ 6t), with fldust(^) = fldust(400) ■ 




The parameter 6t above is the temperature difference relative to the depth of 1730 m (center 
of AMANDA): 

6T(d) = r(c/)- r(1730m). 
The temperature T(K) vs. depth d(m) is parameterized in (jsl as: 

T = 221.5 - 0.00045319 ■ d + 5.822 ■ lO"*" ■ d^. 

The two remaining global parameters, D and E, were defined in [4^1 in a relationship estab- 
lishing a correlation between absorption and scattering: adust(400) ■ 400*^ ~ D ■ bg{4QQ) + E, but 
are not used in this paper 

This work presents the measurement of the values of be{4QQ) and a(400) based on data taken 
at a wavelength of 400 nm and relies on the six-parameter ice model described above to extrap- 
olate scattering and absorption for wavelengths other than 400 nm. 



5. Simulation 

The detector response to ffashing each of the 60 DOMs on string 63 generated a large data 
set that required very fast simulations such that many different sets of the coefficients be(4-0Q) 
and fldust(400) could be compared efficiently with the data. A program called PPC (photon prop- 



agation code, see Appendix A i, was written for this purpose. PPC propagates photons through 
ice described by a selected set of parameter values for be{4QQ) and adust(400) until they reach a 
DOM or are absorbed. When using PPC, no special weighting scheme was employed except that 
the spherical DOMs were scaled up in radius by a factor of 5 to 1 6, depending on the required 
timing precisioij^ and the number of emitted photons was scaled down by a factor of 5^ to 16^, 
corresponding to the increased area of the DOM. 



* Special care was taken to minimize any bias on plioton arrival times by oversizing DOMs. First, we oversize 
DOMs in the direction perpendicular to the photon direction in order to avoid an artificially reduced propagation path 
before reaching the receiver Still, in the worst case, an increase in size by a factor of 16 above to the nominal DOM 
dimensions may introduce a maximum bias of (16 - 1) ■ 16.51 cm / 22 cm/ns=l 1.3 ns towards earlier arrival times (for 
a DOM with radius 16.51 cm and for speed of light in ice of 22 cm/ns). However, on average this error is smaller 
An additional consideration is the overestimated loss of photons if they would get absorbed when entering an oversized 
DOM. Therefore we allow the photons to continue propagating even after hitting an oversized DOM. 
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Figure 7: Angular sensitivity of an IceCube DOM where is the photon arrival angle with respect to the 
PMT axis. The nominal model, based on a lab measurement, is normalized to 1.0 at cos;; = I. The area 
under both curves is the same. 
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Figure 8: Left: optical module acceptance: fraction of photons aiTiving from a direction parallel to the PMT 
axis (at cos;; =1) that are recorded. Note that the acceptance here is meant to include the glass and gel 
transmission and the PMT quantum and collection efficiencies. The acceptance is substantially lower at 
the peak than the roughly 20-25% quantum efficiency of the PMT alone because it is given with respect to 
the photons incident on a cross-section of a DOM, which is larger than that of a PMT. Right: number of 
Cherenkov photons emitted by one meter of the simulated bare muon track (i.e., muon without secondary 
cascades), convolved with the optical module acceptance. The integral under the curve is 2450 photons. 



The relative angular sensitivity of the IceCube DOM was modeled according to the "hole 
ice" description of |9], which is shown in Fig.|7] The "hole ice," a column of ice approximately 
30 cm in radius immediately surrounding the IceCube string, is described by taking into account 
an increased amount of scattering (with effective scattering length of 50 cm) via an empirical 
modification to the effective angular sensitivity curve of the receiving DOM. 

The DOM acceptance is defined as the fraction of photons incident onto the cross-section of a 
DOM that cause a signal in the PMT. This fraction accounts for the losses due to the glass and gel 
transmission, and includes PMT quantum and collection efficiencies. It was calculated according 
to 101 for a DOM of radius 16.51 cm. At 400 nm the DOM acceptance for the photons arriving 
at the PMT along its axis is 13.15%. This corresponds to the nominal angular sensitivity curve 
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of Fig. |2]peaking at 1 .0 for cos t] - 1. Additional considerations, including partial shadowing of 
the DOM surface by the supporting cables, lower this value by 10%. 

The actual number of isolated photoelectrons recorded by a DOM is reduced a further 15% 
because of losses due to the discriminator threshold. The counting efficiency for single photons 
incident on a DOM is thus 13.15% ■ 0.9 ■ (1-0.15) - 10%. The peak of the ampUtude distribution 
for one photoelectron is used to normalize this distribution and is henceforth used as a unit called 
p.e. The discriminator threshold is set at 0.25 p.e. The mean value of the amplitude distribution 
is found at 0.85 p.e., and at ■ 0.85 p.e. for photo-electrons recorded in one sensor. Thus, the 
fraction of charge recorded in multi-photoelectron records is the same as the recorded fraction 
of the number of isolated photoelectrons, 0.85. In a multi-photoelectron dominated situation this 
number can be used to convert from photoelectrons to amplitude in the p.e. unit. The product of 
this value and the two factors listed in the previous paragraph, 13.15% ■ 0.9 ■ 0.85 = 10%, is the 
"effective acceptance," and is applied later (see section[8]l. 

Naturally abundant cosmic ray muons which reach the depth of the detector produce Cherenkov 
light in a broad wavelength spectrum and may be used to test the ice model. For the tests pre- 
sented in section[TT] we simulate the light emitted by muons according to the following method. 
The Cherenkov photons were sampled from a convolution of the wavelength dependence of the 
DOM acceptance with the Cherenkov photon spectrum (see Fig. [8] right) given by the Frank- 
Tamm formula ifioll : 

dN Ina , 
- ■ sm Qr- 



dAdl 

The muon light production is treated via the use of the "effective length" (dl), as described in 
Appendix B The phase refractive index, rip, used in the formula above (defining the Cherenkov 



angle cos 6c = 1 /rip) and the group refractive index, rig, used in calculation of the speed of light 
in the medium, were estimated according to formulas from |l3|l : 

rip = 1.55749- 1.57988- /I -H 3.99993 -/l^- 4.68271 -/l^ -H 2.09354-/ 

rig = rip - (1 + 0.227106 - 0.954648 - A + 1.42568 - A^ - 0.711832 - /). 

The distribution of the photon scattering angle 6 is modeled by a linear combination of two 
functions commonly used to approximate scattering on impurities: 

^(cos 9)^(1- /sl) - HG(cos 6*) + fsi ■ SL(cos 0). 

The first is the Henyey-Greenstein (HG) function 10]: 

HG(cosg)^ ^ ^ 2 ^ with g^icosO), 

2 [1 + g - 2g - cos Oy/^ 

which can be analytically integrated and inverted to yield a value of cos as a function of a 
random number^ uniformly distributed on interval [0; 1]: 



1 

cos 6 - — 
28 



1+.2 



1 2\2\ 



1 +gs 



2-^-1. 



The second is the simplified Liu (SL) scattering function 111 111 : 



l+a /l+cose^" 2g 

SL(cos6») = , with a^—^, g^(cos0}, 

2 \ 2 / l-g 
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which also yields a simple expression for cos as a function of a random number ^ e [0; 1]: 



cos 



= 2-^-1, with 13 



1 +j 



Figure |9]compares these two functions with the prediction of the Mie theory, with dust concen- 
trations and radii distributions taken as described in [4]. The photon arrival time distributions are 
substantially affected by the "shape" parameter /sl (as shown in Fig.fTOb. making it possible to 
determine this parameter from fits to data. 
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Figure 9: (left) Comparison of tlie Mie scattering profiles calculated at several depths of the South Pole ice 
with the Henyey-Greenstein (HG) |4] and simplified Liu (SL) |11] scattering functions. In each, g = 0.943. 

Figure 10: (right) Photon arrival time distributions at a DOM located 125 m away from the flasher, simulated 
for several values of g = (cos 9) and /sl. The difference in peak position simulated with g = 0.8 and g = 0.9 
is of the same order (~ 10 ns) as that between sets simulated with different values of the shape parameter 

/SL. 



A value of ^ = 0.9 was used in this work (cf. g - 0.8 in f?]). A higher value, 0.94, is predicted 
by Mie scattering theory |4, 12], but results in slower simulation and almost unchanged values 
of the effective scattering {be) and absorption {a) coefficients, as shown in [4J. 



6. Likelihood description 

Consider the amount of charge received by DOM / in time bin n when flashing DOM k. 
This charge is measured by taking data with a total photon count of d in rid flasher events and 
a per-event expectation of fi^- This charge is predicted by the simulation with a total photon 
count of s in n, simulated events and a per-event expectation of yUj. Naively one expects the best 
approximations to ji^ and from data and simulated events to be fid - d/tid and //j = s/hj. 

The error in describing data with simulation (i.e., describing fi^ with jj,,) is approximately 
20 - 30% (estimated as described later in section [TOli. One quantifies the amount of disagreement 
between data and simulation in the presence of such an error with ^x]nk- Omitting the indices /, 
n, and k, this is given by: 

2 (In /Id - ln/Zj)2 

^ = 2 • 
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The uncertainty due to this systematic error can be modeled with a probability distribution func- 
tion 

1 -Qnnd-lnfisf 
, exp — r . 

Given that ju^ and yu, are not known, and the measured values are d and s, one formulates the 
likelihood function that describes counts measured in both data and simulation as 

Taking the negative logarithm, this becomes: 

In 5! + /z ^n, - s ln(//,n,) + In <i! + /z^n^ - d InQidnd) + — ^ In^ — + ln( Vlncr) = F. 

2cr^ Us 

The function F(jis,^d) can be easily minimized against jUj and fid, yielding estimates of these 
quantities. To demonstrate this, the derivatives of F are first calculated and set to 0: 



OF i , fid 
sT— ^Hsris-s rln — 



dF , ^ ^ l^d 

l^d^ = IJdnd - d + ^ In — = 0. 

The sum of these Ou^nj + HdUd = s + d) yields an expression of Hd as a function of Hs. Plugging 
it into the first of the above two equations one gets 



dF 1 /irf(jUj) 
<^(Ms,fid(Ms)) = f^sns -s r In 



This equation can be solved with a few iterations of the Newton's root finding method starting 
with a solution to 

s + d 

fls = l^diMs)- Hs=fid = • 

Hs +nd 

At each iteration the value of ju, is adjusted by -///', where the derivative is evaluated as 

/'=nil + ^(— + — ) 

Once the likelihood function is solved for the best values of yU, and yU^, they may be inserted 
into the expression for x^„ ^ above. One can now write the complete function (adding the 
regularization terms Rj described in the next section) as a sum over DOMs i, used in the analysis, 
multiplied by time bins n, when flashing DOMs k: 

^ =L — ^2 — + L ""^^J- 

i,n,k j={r,u} 
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7. Regularization terms 



Two regularization terms (see, e.g., 111311 ') are added to the likelihood function described in 
the previous section. The first term suppresses the fluctuations of scattering and absorption co- 
efficients with depth in under-constrained ice layers. It is formed from terms that are numerical 
expressions for second derivatives of scattering and absorption with respect to the position of the 
ice layer: 



Rr = Yj - 1] - 2 ■ In be[i] + InbeU + i]f 

+(lnfldust[! - 1] - 2 ■ lnadust['] + lnadust[! + 1])^] ■ 



i=2 



Here is the number of ice layers in which b,, and Adust are defined. 

The second term is used to suppress fluctuations in the diagram of Odust vs. b,,, enforcing 
the notion that both are proportional to the dust concentration. It is constructed as an excess 
of the sum of distances between the consecutive points (In bg. In fldust) over the shortest distance 
connecting the end points: 

N-l 



R„^-D{l,N) + YjDijJ+l), 



where £>(;'i, j2) = -J(inbe[ji] -^be[j2])^ + (InadustU'i] - InadustU'a])^- 

The points (In b^, In fldust) are sorted by the value of In b^, + In fldust and shown in the above sum 
with the index j[i]. 

Both terms affect the resulting scattering and absorption coefficients by on average less than 
2% at detector depths at their chosen strengths a,-„. Deviations larger than this, up to 19% were 
observed in the region of particularly dusty ice around the depth of 2000 m. The size of the effect 
has been verified by re-running the fits without including the terms. The regularization terms are 
likely to become more important if the thickness of ice layers (10 m in this work) were chosen 
to be much smaller than the spacing between DOMs on a string (17 m). 



8. Fitting the flasher data 

The six horizontal LEDs within a single DOM flashing at maximum brightness and width 
nominally emit about 4.5 ■ 10'" photons 1 1]. After accounting for the effective DOM acceptance 
(as explained in section |5]), these photons result in a charge amplitude of 4.5 ■ 10^ p.e., which 
henceforth is traced as 4.5 ■ 10'' "photons" that each result in an amplitude of I p.e. Using a 
DOM size scahng factor of 16, only 1.76 ■ 10^ photons need to be simulated (16^ = 256 times 
fewer). 

Simulating 9765625 photons, with a scaling factor of 16, corresponds to 2.5-10' photons sim- 
ulated for actual-size receiving DOMs, or 2.5 ■ 10'" real photons leaving the flasher DOM (after 
accounting for the effective acceptance of the receiving DOM). This is defined as a "unit bunch" 



of photons, which is simulated in approximately 1 second on a single GPU (see Appendix Ai. 

In the following discussion, a "photon yield factor" (py) is the number of unit bunches that 
corresponds to a given number of real photons. For instance, 4.5 ■ 10'" photons emitted by a 
flasher board correspond to a photon yield factor of py = 1 .8. 
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Data from all pairs of emitter-receiver DOMs (located in the same or different ice layers, 
amounting to about 38700 pairs) contributed to the fit to 200 ice parameters (scattering and 
absorption in 10 m layers at detector depths of 1450 to 2450 m). Two functions were used in 
fitting the data: xl was constructed with one term from each emitter-receiver pair (using the total 
recorded charge) andxj was constructed with the recorded charge split in 25 ns bins. Although 
X^ more completely used the available information, x^ was found to be somewhat more robust 
with respect to statistical fluctuations in the simulated sets and was faster to compute. Thus, x^ 
was used in an initial search for a solution, with;^'^ applied in the final fits. 

Both be{4QQ) and adust(400) are roughly proportional to the concentration of dust (this would 
be precise if the dust composition in the ice were the same at all depths). This motivates the fol- 
lowing simplification in the initial search for the minimum ofx^' in each layer both bg(4QQ) and 
adust(400) are scaled up or down by the same relative amount, ranging from 1-40%, preserving 
their ratio to each other 



Starting with some initial values of be(400) oc Odusti^OO) and some /?,., foff, /sl^ 
Using xi find best values of beiAQQi) ~ adust(400) 
Using xj find best values of py, fotf, /sl, ^sca, a-^bs- 

Py'. photon yield factor 

foff: global time offset for all flasher pulses 

fsL- shape parameter of the scattering function 

asca^ scaling of scattering coefficient table 

ffabii: scaling of absorption coefficient table 
repeat this box until converged (~ 3 iterations) 
Using x^, refine the fit with be{4-Q0) and fldust(400) fully independent from each other. 

Table 1 : Flow chart of the global fit procedure to ice/flasher parameters. 

Starting with homogeneous ice described by be(4QQ) = 0.042 m"' and fldust(400) - 8.0 km"' 
(average of 10] at detector depths), the minimum of x^ is found in about 20 steps. At each 
iteration, the values of ^7^,(400) and fldust(400) are varied across consecutive ice layers, one layer 
at a time. Five flashing DOMs closest to the layer in which the properties are varied are used to 
estimate the variation of the;^'^. FigurefTTIshows fitted ice properties after each of 20 steps of the 
minimizer 

The search for the minimum of x^ is performed next in the parameter space of the overall 
time off'set from the flasher start time (foff), photon yield factor (py), shape parameter (/sl) of 
the scattering function (see sectionjSj, and scaling coefficients ajca and o-abs applied to the depth 
tables of /7,(400) and adust(400). 

The fee(400) and fldust(400) of the solution are scaled with coefficients Csca and o-abs to pro- 
duce the likelihood profiles shown in Fig. [12] From this figure, it is apparent that using the timing 
information is necessary to resolve both ^7^,(400) and fldust(400). The minimum of has an elon- 
gated shape, and the direction of its longest extension is determined. The point on the line drawn 
in this direction through the minimum of Xq is chosen to minimize the xj- The global scaling 
factors ffsca and Oabs corresponding to this point are used to rescale the starting "homogeneous 
ice" values of ^7^(400) and fldust(400). The entire procedure is then repeated. 

The solution is finally refined by varying be(4QQ) and fldust(400) at each step of the xj mini- 
mizer four times (combinations of and adusti^^dust, with6be/be and ^adust/^dust - 1-2%). 
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depth [ m ] iteration of the minimizer 

Figure 11: Left: values of /^^ (400) and a(400) vs. depth after 20 steps of the minimizer. The black curve 
shows fitted values after the last step of the minimizer. Right: global x^q values achieved after each step of 
the minimizer. The starting "homogeneous ice" value is 1.34- 10^. Regularization tenns R,.„ use the scale on 
the right. Also shown are the Poisson terms for simulation and data (llhsj) and the full likelihood including 
all terms (llhtot). The;j'^ changes suddenly when the number of simulated flasher events is increased, but 
ultimately decreases as the minimizer steps through the iterations. Note that for iteration steps 1-10, only 
1 flasher event is simulated. For steps 11-15 and steps 16-20, 4 events and 10 events, respectively, are 
simulated. 
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Figure 12: Likelihood functions in the vicinity of their minima constructed using only charge information 
(left), and using full timing information (right). The values are shown on a log scale (with colors and 
contours). The ranges of values shown are: x\ = 1-43 ■ 10'' to 1.51 ■ 10^ (left) and = 1.05 ■ lO' to 
4.01 ■ 10^ (right). 



The entire procedure described above is also outlined in Table[T] 

The best fit is achieved for py - 2.40 ± 0.07. Since the best value of py is also calculated 
by the method above, the resulting table of be{4-QQ) and aciust(400) values is independent of a 
possible constant scaling factor in the charge estimate or the absolute sensitivity of a DOM. The 
best fit values of the other parameters are foff = 1 3 + 2 ns and fsi = 0.45 + 0.05 (see Fig.flJ]). The 
typical agreement of data and simulation based on these parameters is demonstrated in Fig.|2]5- 
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Figure 13: Behavior of and^^ in the vicinity of the fit minimum vs. r^ff, p,,, and /sl. All plots are shown 
on a linear scale. Horizontal dotted lines show the ±1(T range due to purely statistical fluctuations in the 
simulation estimated for the best-fit model. The minimum in /off and /sL is a feature of but not Xg ■ 



9. Dust logger data 

Several dust loggers fs*] were used during the deployment of six of the IceCube strings result- 
ing in a survey of the structure of ice dust layers with an effective resolution of approximately 
2 mm. These layers were then matched across the detector to provide a tilt map of the South 
Pole ice, as well as a high-detail average dust log (a record of a quantity proportional to the dust 
concentration vs. depth). Additionally, the East Dronning Maud Land (EDML, see [SJ) ice core 
data were used to extend the dust record to below the lowest dust-logger-acquired point (taken at 
a depth of 2478 m). 

The table of dust layer elevations (the tilt map) taken from |5] provides the layer shift (relief) 
from its position at the location of a reference string, at a point distance r away from this string, 
along the average gradient direction of 225 degrees SW (see Fig.[T4li. The z-coordinate of a given 
layer at r is given by Zr - Zq + relief(zo, r). Between the tabulated points, z,- was calculated by 
linear interpolation in zo and r. The equation was solved by simple iteration resulting in a table 
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of zoizr) - Zr VS. Zr giveii at several points along the gradient direction. Combined with the dust 
depth record at the location of the reference stting (at r = 0), this yields a complete description of 
the dust profile in and around the detector (assuming that the concentration of dust is maintained 
along the layers). 
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Figure 14: Top: extension of ice layers along the average gradient direction. Relief is amplified by a factor of 
3 to enhance the clarity of the layer structure. The lowest layer shown exhibits a shift of 56 meters between 
its shallowest and deepest points (which is the largest shift of all layers shown in the figure). Bottom: 
comparison of the average dust log with the effective scattering coetficient bei400) measured with the flasher 
data. 



The correlation between the effective scattering coefficient measured with the IceCube flasher 
data and the average dust log (scaled to the location of string 63) is excellent, as shown in Fig. 
[T4l All major features agree, with well-matched rising and falling behavior, and are of the same 
magnitude. Some minor features are washed out in the flasher measurement. 

With an established correlation to the average dust log, the EDML-extended version of the 
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log was utilized in constructing an initial approximation (replacing the "homogeneous ice") used 
by the fitting algorithm described in section |8] This resulted in the recovery and enhancement 
of several features in the scattering and absorption vs. depth that were previously unresolved. 
Additionally, the solution is now biased towards the scaled values of the extended log (instead 
of the somewhat arbitrary values of the initial homogeneous ice approximation) in the regions 
where the flasher fitting method has no resolving power, i.e., above and below the detector 



10. Uncertainties of the measurement and final result 



To study the precision of the reconstruction method, a set of flasher data was simulated with 
PPC (250 events for each of the 60 flashing DOMs). The agreement between the simulated 
and reconstructed ice properties is within 5% at depths of the instrumented part of the detector 
(see Fig. [TsT l. Due to the dramatically lower number of recorded photons in the layer of ice 
containing the most dust, at around 2000 m (the dust peak), additional simulation was necessary 
to reconstruct the local ice properties: 250 events per flasher were used within the dust peak, 
whereas only 10 events per flasher were used elsewhere. Up to 250 simulated events per flasher 
were used to achieve the best possible precision of the final result shown in Fig. [16] 
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Figure 15: Reconstructed ice properties in black for simulated flasher events with input ice properties in red. 
The ice properties in the dust peak are reconstructed correctly with 250 simulated events per flasher. The 
blue dashed curve shows the result achieved with only 10 simulated events per flasher. 



This verification approach was used to quantify the uncertainty in the measured values of 
bi,(4QQ) and fl(400) due to the lack of knowledge of the precise flasher output timing profile. 
Reconstructing the simulation, which used the 62 ns-wide rectangular shape of the flasher pulse, 
with a hypothesis that all photons are emitted simultaneously at flasher start time, leads to max- 
imal systematic shifts in the obtained effective scattering and absorption coefficients of roughly 
6.5%. 

Several pulse extraction methods, with and without correcting for PMT saturation (using the 
saturation model of ||3l), were tested for extracting photon hits from the flasher data, and the ice 
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properties were reconstructed for each and compared. This provided an estimate of about 4% for 
the uncertainty in the measured ice properties (effective scattering and absorption coefficients) 
due to detector caUbration and pulse extraction (in waveforms of up to 1000 photoelectrons). 
We note that reconstructing the data with the azimuthally symmetric vs. 6-fold "star" pattern of 
flasher LED light leads to no discernible difference in the resultant ice properties. Further, since 
the DOMs on the flashing string are not used in the fits, the difference between the ice properties 
reconstructed for nominal or hole ice angular sensitivity models is negligible. 

Finally, the uncertainty due to statistical fluctuations in the sets simulated during the recon- 
struction procedure are estimated at roughly 5-7%. This uncertainty could be reduced with more 
simulated events per flasher (at least 10 were simulated for each configuration, compared to 250 
events present in data). However, given that the entire fitting procedure currently exceeds 10 days 
of calculation to produce a result, the number of simulated events represents a Umiting constraint. 

The effective scattering and absorption parameters of ice measured in this work are shown 
in Fig. [16] with the +10% grey band corresponding to ±lcr combined statistical and systematic 
uncertainty at most depths (values of 6.5%, 4%, and 5 -7% explained above are added in quadra- 
ture to result in a total uncertainty estimate of 9.7+0.6%). The uncertainties may be somewhat 
larger than this average value in the layers of dirtier ice, since many of the detected photons are 
likely to spend more of their travel time in the adjacent layers of cleaner ice (thus resulting in a 
weaker constraint of the properties of a dirtier ice layer). The uncertainty increases beyond the 
shown band at depths outside the detector, above 1450 m and below 2450 m. 

Figure [16] also shows the AHA (Additionally Heterogeneous Absorption) model, based on 
the ice description of |4], extrapolated to cover the range of depths of IceCube and updated with 
a procedure enhancing the depth structure of the ice layers. The AHA model provided the ice 
description used for Monte Carlo simulations of IceCube data prior to this work. 

How well we fit the ice properties is limited by our ability to properly simulate all propagation 
and instrumental effects. Not all effects are accounted for, as it appears, in the analysis presented 
here, despite the simplicity of the physics model involved. In order to estimate the error in the 
description of the data with the fitted model, we created a histogram (see Fig.fTTTi of the ratio of 
simulation to data for sufficiently large charges, minimizing statistical effects. The width of this 
histogram, estimated to be around 30% of the received charge, represents the "model error" and 
enters the fit as a parameter in the likelihood function (see section|6]l. 

It is not possible to translate all of this model error into the uncertainties in the measured 
parameters since we can only estimate uncertainties from the known causes (e.g., by varying the 
parameters of the PMT saturation model). During our investigation of the discrepancy demon- 
strated in Fig. [17] we found evidence of different propagation properties of photons in different 
directions inside the detector. Nevertheless, the uncertainties given above are still applicable to 
average (over all directions) values of effective scattering and absorption. The resulting situa- 
tion compels us to report both the parameter uncertainties (~10%), and the average model error 
(~30%, when describing the flasher data as shown in Fig. [17). 

11. Comparison of full-detector simulation to muon data 

To investigate the accuracy of the resultant ice model in describing actual IceCube data, anal- 
yses were performed that compared experimental muon events with simulation. Atmospheric 
muons are a source of physics events for IceCube but represent a background for neutrino anal- 
yses. In the 2008 40-string configuration, atmospheric muons triggered IceCube at a typical 
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Figure 16: Values of the elfective scattering coefficient ^^(400) and absorption coefficient a(400) vs. depth 
for a converged solution are shown with a solid line. The range of values allowed by estimated uncertainties 
is indicated with a grey band around this line. The updated model of 1 4] (AHA) is shown with a dashed 
line. The uncertainties of the AHA model at the AMANDA depths of 1730 ± 225 m are roughly 5% in 
and roughly 14% in a. The scale and numbers to the right of each plot indicate the corresponding effective 
scattering Ijhg and absorption 1/a lengths in [m]. 



rate of 1 kHz, and therefore a large statistical data set was available for comparisons between 
measured muon data and simulations of cosmic ray induced muons. The simulations are based 
on the assumed propagation of optical Cherenkov photons through the ice but also depend on 
assumptions that include the energy, multiplicity, and angular distribution of the muons. 

The simulation chain begins with the production of atmospheric muons from cosmic ray air 
showers using the CORSIKA software lll411 . followed by propagation of the muons with muon 
Monte Carlo (MMC) 1150 and generation of photons according to a Cherenkov spectrum and 
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Figure 17: Histogram of the ratio of tlie average charge per emitter-receiver pair in simulation to data. The 
lower histogram has one entry per pair, based on the averaged total charge. The upper one has one entry per 
time bin. Pairs and time bins are used only if the average charge in the data is more than 10 photoelectrons. 
The widths of the fitted log-normal distributions are 0.29 and 0.31, respectively. 



their propagation with photonics |9] or PPC lll6ll . Finally the photons are detected and digitized 
by the DOM simulator To compare different ice models and photon propagation techniques, 
only the parameters relevant for the photon propagation are varied in simulation, while all other 
settings remain fixed. 



11.1. Delta T distributions 

A relatively generic method to compare ice models and examine specific ice properties de- 
scribed here utilizes DeltaT distributions. DeltaT is defined as the time difference between first 
hits on adjacent DOMs on the same string. A positive DeltaT represents a photon that strikes the 
upper DOM followed by a photon strike of the DOM directly below. This method permits close 
investigation of basic photon timing information without requiring ice-model-dependent muon 
reconstruction techniques. The distribution of DeltaT values for downgoing muon data taken 
with the 59-string detector configuration during September 2009 is shown in Fig. [18] The tails of 
this distribution consist of relatively long-lived photons and contain information about the bulk 
ice properties, such as scattering and absorption. On the other hand, the peak of the distribution 
consists of photons that travel from source to DOM with few scatters (i.e., "direct" photons), 
and is relatively invariant to the depth-dependent bulk ice properties. Figure [19] illustrates this 
relationship throughout all detector depths. 

Full-detector Monte Carlo simulation was generated for different ice parameters to examine 
their effects on the shape and height of the peak in the DeltaT distribution. Figures l20li23] show 
the peak shape for data and various simulation models. The description of the ice denoted as 
SPICE2x was an intermediate model in this analysis, and is characterized by similar scattering 
and absorption lengths to those of the SPICE Mie model, which is the final result. In SPICE2x, 
a Henyey-Greenstein (HG) scattering function is used instead of a linear combination of the HG 
and SL functions. Additionally, SPICE2x has an average scattering angle of g - (cos 6) - 0.8 
instead of 0.9 (used in the final result), and lacks the global flasher time offset parameter used 
in the fit of SPICE Mie. In all of the permutations of the ice properties examined, the only 
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Figure 18: Left: DeltaT distribution for muon data. Tlie cutoff at ± 1000 ns is due to the coincidence trigger window 
where data from a triggered DOM will only be read out if an adjacent or next-to-adjacent DOM also triggers within a 
time window of 1000 ns. Right: A zoom of the peak of the distribution. The peak is shifted towards positive times 
because it is dominated by direct photons from downgoing muons, which are detected first by the upper DOM and then 
the lower DOM. The shift roughly corresponds to the muon flight time between DOMs. 
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Figure 19: DeltaT distributions for DOMs binned in 20 m depths. Widths of the peaks and tails were extracted and plotted 
vs. depth for the entire detector. Left: The full width at 5% of the maximum, corresponding to the width of the tails, 
shows a strong depth dependence similar to the dust logger data and the derived scattering and absorption parameters. 
Right: The full width at half maximum (FWHM) shows very little depth dependence. The FWHM was computed by 
multiplying the number of bins and the bin width, resulting in the discrete level structure in the plot. 



parameters that significantly clianged tlie sliape of tlie peak were the hole ice scattering, scattering 
function composition, and the time offset parameter. 



77.2. Event geometry and time residuals 

The simulation data for this study were produced for the IceCube detector in its 40-string 
configuration and is compared to data taken in August 2008. This corresponds to roughly 10% 
of the yearly experimental data of IceCube. 

For a generic comparison, it is preferable to use an unbiased data sample. For such purposes, 
IceCube operates a Minimum Bias filter stream that selects every event that was recorded by the 
DAQ independently of the satisfied trigger condition. A prescale factor of 2000 was applied to 
data events to comply with bandwidth requirements before sending data north via satellite. This 
analysis used events that passed a DOM multipUcity condition of at least 8 DOMs within a 5 //s 
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Figure 20: The peak region of the DeltaT distribution for the SPICE2x model shows a lack of direct photon hits compared 
to the data. Neither increasing the amount of forwai'd scattering (by setting g = 0.95) nor increasing the bulk ice scattering 
by 20% significantly changes the peak height or shape. 



time window that register a hit in coincidence with one of their vertical neighbors (within 1 fus). 
From this data stream, events that had sufficient recorded information to perform reconstructions 
of reasonable accuracy were selected. The selection criteria were based on the zenith angle 
{6 < 90°) and the likelihood minimum of the standard angular fit lll7ll divided by the degrees 
of freedom (reduced log-likelihood, rlogl < 8). The resulting median angular resolution of this 
event sample was better than 2° with a passing rate of roughly 15% of the initially recorded 
data. The comparisons shown in Figs. I24ll26l are based on 130 million events. The absolute 
normalization between experimental data and simulation was affected by large uncertainties, but 
for the purpose of this analysis all distributions were normalized to unity, and the focus was 
placed on differences in shape between the curves. 

A basic test to examine the influence of ice parameters on the simulation is to compare depth- 
dependent variables since the optical ice properties characteristically vary with depth. Figure l24l 
shows the distribution of hit DOMs. The peaks correspond to regions with clearer ice, and the 
valleys indicate ice containing more dust. The simulation that uses a combination of the SPICE 
Mie ice model and the PPC photon propagator shows a significant improvement in agreement 
with the experimental data. The ratio between Monte Carlo and experimental data histograms 
is almost flat, except for a dip around DOM 36, which is a region of high dust concentrations 
and therefore poor statistics. Figure |25] shows the distribution of the z-coordinate of the center 
of gravity of all hit DOMs for each event. A marked improvement in the agreement between 
experiment and simulations, in particular in the deep ice, is observed for the SPICE Mie and 
PPC Monte Carlo. 

Similar to the timing distributions for the flasher events that were used to extract the ice 
properties, the distribution of arrival times for Cherenkov light from muons can be investigated. 
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Figure 21: The peak region of the DeltaT distribution shows sensitivity to the hole ice description. The hole ice is 
modeled as a vertical column of ice with a higher concentration of air bubbles, resulting in a local scattering length of 50 
cm. Simulations that assume no contribution to scattering due to hole ice and those with three times the nominal bubble 
concentration in the hole ice are shown. The hole ice is thought to increase the number of direct photon hits because the 
increased scattering in the region of the hole ice causes more photons from downgoing muons to be locally up-scattered 
into the PMT (which faces downward in the DOM), effectively altering the angular sensitivity of the DOM. 



Here we study the distribution of time residuals, which are calculated by subtracting the expected 
geometrical arrival time for unscattered Cherenkov photons (based on the estimate of the recon- 
structed track) from the actual arrival time of the detected photons. If the track is reconstructed 
accurately, the time residual is caused by scattering. The slope of the time residual distribution 
is strongly correlated with the optical absorption length. The new simulation shows an improved 
overall description of the experimentally observed timing residuals, see Fig.|26] The plot of the 
ratio between Monte Carlo and experimental data histograms shows an almost flat behavior for 
the most relevant time interval up to 1 fis. Note that the distribution is summed over depth so 
discrepancies at specific depths may remain. 



12. Conclusions and outlook 

The precise modeling of the optical properties of the South Pole ice is crucial to the analysis 
of IceCube data. The scattering and absorption coefficients of ice (averaged in 10 m depth bins) 
were obtained from a fit to the in-situ light source data collected in 2008. The result is shown 
in Fig. [16] and also presented in [Appendix C| The sum of the statistical and systematic uncer- 
tainties in the measured values of the effective scattering and absorption coefficients inside the 
instrumented volume of ffie IceCube detector was estimated to be less ffian 10%. 

r— t 

This analysis builds on the concepts developed in [4], and relies on the wavelength depen- 
dence determined there. Unlike in |4], this analysis uses a global fit to simultaneously describe 
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Figure 22: The peak region of tlie DeltaT distribution for SPICE Mie, comparing the full model (fst = 0.45) to the 
model with only the HG scattering function (i.e., setting fst = 0). The obsereed effect is thought to be caused by the 
higher probability of photons scattering through intermediate angles of 20° -40° . Even though the typical muon-to-DOM 
distance is small compared to the effective scattering length, photons are more likely to scatter at larger angles and 
therefore to be detected. 




Figure 23: The peak region of the DeltaT distribution comparing the final SPICE Mie fit result to the previous AHA 
model and the muon data. The fit to the data is significantly improved with the SPICE Mie model. 
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Figure 24: Distribution of hit DOMs. Tlie label "OM_N" is the number of the DOM on the stiing, ranging 
from 1 at the top of the detector to 60 at the bottom. The curves are normalized to one event for a better 
comparison of the shape. The plot on the right shows the ratio between simulation and data. 




600 



Figure 25: Distribution of the z-coordinate of the center of gravity of all hit coordinates for each event. The 
curves are normalized to one event for a better comparison of the shape. The plot on the right shows the ratio 
between simulation and data. 



all of the light source data. It also uses significantly more data than [4], both in terms of the 
number of registered photons and the number of emitter-receiver pairs. 

The high quality of the fit was ensured by careful selection of the likelihood function that 
quantified the differences between data and simulation within a single model of statistical and 
systematic uncertainties. In the course of the investigation we found that determining the shape 
of the scattering function and incorporating the ice layer tilt was necessary to achieve better 
model agreement with the data. 

We are aware of some systematics issues that are not yet entirely understood (and will be the 
subject of further studies). One notable omission from this work is the direct simulation of the 
photon propagation in the columns of ice refrozen around the IceCube strings. A study of the 
slight anisotropy hinted at in this report (section [10) is the subject of a forthcoming publication. 
Additionally, we have not yet analyzed the data from the LEDs with wavelengths other than 
400 nm (which were installed during the final IceCube deployment season in 2010). Thus, the 
new ice model presented here relies on the previously-established wavelength dependence of 
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Tres [ns] 



Figure 26: Distribution of tlie time residuals: time delay due to scattering of photons aniving from the recon- 
structed muon track in data and simulation. The curves are normahzed to one event for a better comparison 
of the shape. The plot on the right shows the ratio between simulation and data. 



effective scattering and absorption coefficients. However, we are encouraged by the significantly 
improved agreement between data and simulation when using the new ice mode obtained in this 
analysis. 



Appendix A. Photon Propagation Code 



Four different versions of the program (available from 111611 ) were written: one in C++, an- 
other in assembly (for the 32-bit i686 with SSE2 architecture), and a version that employs the 
NVIDIA GPUs (graphics processing units) via the CUDA programming interface |18]. The 
fourth version uses OpenCL lUsll . supporting both NVIDIA and AMD GPUs, and also multi- 
CPU environments. The relative performance of these different implementations (for simulating 



both flashers and Cherenkov light from muons) is compared in Table IAT2 



C++ 


Assembly 


CUDA 


OpenCL 


flasher 1.00 


1.25 


147 


105 


muon 1 .00 


1.37 


157 


122 



Table A. 2: Speedup factor of different implementations of PPC compared to the C++ version. C++ and Assembly code 
was run on 1 core of Intel 17 920 (2.66 GHz) CPU. CUDA and OpenCL code was run on 1 GPU of NVidia GTX 295 
video card. 



The writing of the GPU version of PPC was prompted by a similar project called i3mcml ifioll . 
which showed that acceleration by factor of 100 or more compared to the CPU-only version was 
possible. After demonstrating impeccable agreement between test simulation sets made with the 
C++, assembly, and GPU implementations of PPC, and with i3mcml, the CUDA version of PPC 
was chosen for the analysis of this work on a GPU-enabled computer with i7 920 (2.67 GHz) 
CPU and 3 GTX 295 NVIDIA cards (6 GPUs). 
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Appendix B. Muon and cascade light production 

The light yield of the muon and all of its secondaries (ionization losses and delta electrons, 
bremsstrahlung, electron pair production, and photonuclear interaction 115]) with energies below 
500 MeV is parameterized in [20l!|by substituting the length dl of the Cherenkov-light-emitting 
segment of a bare muon of energy E with the "effective length" 

dle« = dl ■ (1.172 + 0.0324 ■ log,(£ [GeV])) . 

The parameterization given above was used in the muon studies of section (TT] However, we are 
aware of an updated parameterization of [21] and list it here for completeness: 

t/4ff = dl ■ (1.188 + 0.0206 ■ log^(£ [GeV])) . 

The light yield of cascades is also parameterized in lEoll via use of the "effective length": 

dkff = 0.894 ■ 4.889/p m/GeV ■ E [GeV] for electromagnetic cascades 
dies = 0.860 ■ 4.076/p m/GeV ■ E [GeV] for hadronic cascades. 

These formulae were derived for muons in water, but are given here for propagation in ice (p - 
0.9216 is the ratio of the densities of Iceland water). This work relies on newer parameterization 
of the cascade light yield of 1 22 

dlgfi =5.21 m/GeV ■ 0.924/p ■ E [GeV] for electromagnetic cascades 
dkff = F ■ 5 .2 1 m/GeV ■ 0.924/p ■ E [GeV] for hadronic cascades. 

Here F is a ratio of the effective track length of the hadronic to electromagnetic cascades of the 
same energy E. It is approximated with a gaussian distribution with the mean (F) and width cr/r: 

<F)= 1 -(£ [GeV]/£o)""' ■ (1 -/o), £o = 0.399, m = 0.130, /o = 0.467, 
o-F = (F) ■ 6o ■ logio(£ [GeV])-^, 6o = 0.379, y = 1.160. 

In order to properly account for the longitudinal development of the cascade, the distance 
from the start of the cascade to the point of photon emission is sampled from the following 
distribution (ignoring the LPM elongation) 1I20II : 

I = Lrad ■ r{a)/b, Lrad = 35.8 cm/p, 

where r(fl) is a gamma distribution with the shape parameter a. Parameters a and b are given by: 

a - 2.03 + 0.604 ■ logg(£'), b - 0.633 for electromagnetic cascades 

a = 1 .49 + 0.359 ■ log^X^), b = 0.772 for hadronic cascades. 

All photons are emitted strictly at the Cherenkov angle with respect to the emitting track 
segment. These, except for the bare muon itself, are assumed to be distributed according to 

dl/dx ~ exp(-b ■ x") ■ x""', with x - 1 - cos(O). 

The coefficients a - 0.39 and b = 2.61 were fitted to a distribution of 100 GeV electron cascades 
from lEoll (see Fig. lB.27] i and are fairly constant with energy and are used to describe the hadronic 
cascades as well. 



'The formula 7.97 contains a typo; however, the caption within Fig. 7.56 (B) is correct, with LOG(E) understood as 
In(£) EE log,(£). 

^Taken at the center of IceCube (depth of 1950 m, temperature -30.4° C); cf. p = 0.9167 at 1 atm. at 0° C. 
'The axis labels in Fig. 3.2 are correct; formula 3.4 needs to be corrected as in this text. 
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Figure B.27: Fit to the angular track-length distribution for 100 GeV electron cascades. The simulation line 
in black is taken from figure 7.44 of 120] , the fit in green is to the function given in this text. 



Appendix C. Table of results 
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Table C.3: Effective scattering length i/bg and absorption length i/a at 400 nm vs. depth given at the x,y coordinates 
coiTesponding to the center of IceCube array. This, together with the value of the shape parameter of the scattering 
function, fsi = 0.45 constitues the "SPICE Mie" model. Additional parameters that this model depends on that were 
derived elsewhere are parameters a, k. A, and B of the six-parameter ice model 3, and ice tilt map of (sj. 
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